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I . INTRODUCTION 



Guided missiles can be classified into four categories 
depending on launch and target position characteristics. The 
categories are Air-to-Air, Air-to-Surface, Surface-to-Air , and 
Surface-to-Surface missiles. 

Another classification among missiles is the guidance 
system of the missile. The missile can be command or homing 
guidance. 

In the command guidance system the missile and target are 
continuously tracked and guided from one or more friendly 
vantage points, and the necessary path for intercept is 
computed and relayed to the missile. 

In the homing guidance system, the missile has a homing 
device onboard which can detect the target and gives the 
necessary path directions for intercept to the missile. The 
homing missile is further subdivided into classes having 
active, semiactive, and passive guidance systems. Active 
detection is when the missile illuminates the target, i.e., 
with a radar, and receives the reflected signals. Semiactive 
detection is when the target is illuminated by a source other 
than the missile and the missile receives the reflected 
signals. Passive detection is used when the target is the 
source of energy, and the missile detects signals that 
propagate from the target. 
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Each of the missiles in the above categories will employ 
one or more of the three guidance laws. These laws are Pursuit 
Guidance, Line-of-Sight Guidance, and Proportional Guidance. 
The first portion of the missile flight path may use one of 
the guidance laws but the terminal phase of flight may be best 
suited for another. 

The present work addresses the design and evaluation of a 
semiactive Surface-to-Air missile using Proportional 
Navigation as the guidance law. A ground based target tracker 
will also developed with the target deviations in position and 
velocity relayed to the missile. 

Chapter II presents a description and comparison of the 
three different guidance laws. The Proportional Navigation 
guidance law will also be developed. In Chapter III the 
missile and target flight path models will be developed using 
the concepts of Chapter II and computer simulation studies 
will be performed. Chapter IV consists of a Luenberger 
observer design, and an evaluation of the estimator, and the 
guidance law over a range of conditions will be conducted. 
Chapter V consists of the development of the ground target 
tracker using the theory of the Kalman Filter and again, 
computer simulations will be included to determine the 
accuracy of the target tracker. 

All computer simulations are developed and conducted using 
the Matrix Laboratory (MATLAB) language. 
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II. MISSILE GUIDANCE 



A. GUIDANCE LAW SELECTION 

The selection of a guidance law is a pre-requisite for 
determining the initial calculations for the model. The 
missile guidance system measures the error between the 
missile's actual and desired course, computes the corrections 
necessary to reduce the error based on the guidance law 
selected, and gives commands to the autopilot to activate the 
controls required to achieve acceptable intercept of the 
target. The miss distance and the acceleration required by the 
missile are functions of the guidance law. 

1. Pursuit Guidance 

The pursuit guidance law is illustrated in Figure 1 
and is described as having the missile velocity vector 
directed toward the target at all times. The missile is always 
heading along the line-of-sight from the missile to the 
target. This guidance law is effective against slow moving 
targets, but the missile may lack sufficient maneuverability 
against fast moving maneuverable targets. 

2 . Line-of-Sight Guidance 

Line-of-sight guidance is used in a beam-rider type 
missile and is illustrated in Figure 2. This guidance law 
requires that the missile remain on a line joining the target 
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TARGET 1234 




Figure 1. Pursuit Guidance Trajectory 

and the target tracker. The purpose of the target tracker is 
to maintain the antenna boresight pointing at the center of 
the reflecting area of the target. This guidance scheme 
normally requires a dedicated fire control system from launch 
to intercept [Ref.l], 

3. Proportional Navigation Guidance 

This guidance law requires that the missile travel in 
such a way that its own rate of turn is proportional to the 
rate of turn of the line-of-sight from the missile to the 
target. Figure 3 illustrates the proportional guidance scheme 
in which the rate of change of the missile heading is made 
proportional to the rate of change of the line-of-sight 
between the missile and the target. The fixed or variable 
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Figure 2. Line-of-Sight Trajectory 

multiple between the missile rate of turn and the rate of turn 
of the line-of-sight is called the navigation ratio (NR) . The 
proportional navigation guidance law attempts to generate an 
acceleration command perpendicular to the line-of-sight. One 
way to achieve this could be lateral acceleration coupled with 
angular or angular rate commands to place the acceleration 
perpendicular to the line-of-sight. The advantage of this 
guidance law is in its effectiveness against maneuvering 
targets. Since proportional navigation guidance anticipates 

i 

the targets future position, it can attain a higher degree of 
responsiveness over the other two guidance laws. 

Figure 4 illustrates the proportional navigation 
scheme [Ref. 2]. If the seeker head of the missile follows the 
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Figure 3. Proportional Guidance Trajectory 

A. = A. = - 4 - (u R) + w R (2.1) 

T * dt 

target, the target acceleration perpendicular to the line-of- 
sight will equal the acceleration of the R vector, where 
R = missile-target line-of-sight vector 

R = closing rate along R 

u> = angular rate of change of R 
A t = target acceleration perpendicular to R 
A r = tangential acceleration of vector R 
The term (uR) represents the vectorial acceleration of 
R and the term d/dw(wR) represents the rate of change of the 
tangential velocity. A missile acceleration, A M , equal to the 
target acceleration Aj-, at this point will make the line-of- 
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sight parallel to its original direction. Since the velocity 
R ' is along the vector R, a missile/target intercept is 
assured. Therefore, from Equation (2.1), the executed missile 
acceleration commands should be 

A m = 2 (ufl) +((j>R) (2.2) 

Since the direction of the velocity vector cannot be 
directly controlled, proportional navigation is achieved by 
'controlling the commanded missile acceleration (acorn) . 

acom =V *Y (2.3) 

where 
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V M = the missile velocity 
7 = the rate of change of the velocity vector 
Implementation of proportional navigation results in 
the following guidance law 

acom - NR V H 6 (2.4) 



where 

d = the rate of change of the line-of-sight 
With this definition, the equation for the rate of 
change of the velocity vector can be written as 

y = NR d (2.5) 

B. PROPORTIONAL NAVIGATION KINEMATICS 

From Figure 5 the following equations of motion are 
obtained. The three dimensional linear model will be developed 
in Chapter III, but for simplicity and ease of understanding, 
the fundamental equations will be first developed in the x and 
y planes. The velocity vector V M is at an angle y M from the 
established reference line. From the geometry of the problem, 
the missile flight path angle can be easily determined. 
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Figure 5. Intercept Geometry for Proportional Navigation 

y M = arctan (2.6) 

where Vmx and Vmy are the components of the missile velocity 
vector in the x and y directions, respectively. The target 
flight path angle is found from the same geometry and is 
expressed as 

1 T - arctan (2 . 7 ) 

( 

F 

where and Vjy are the target velocities in the x and y 
directions, respectively. 

The missile-target line-of-sight vector is R. 
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R = [ ( X T - X H ) 2 * ( Y t - Y m ) 2 ] 1/2 (2.8) 

where X T and Y x are the x and y coordinates of the target and 
X M and Y m are the x and y coordinates of the missile. 

The magnitude of the velocity vector for the missile and 
target (V T ) can be expressed as 

= [ ( Vmx ) 2 + ( V m ) 2 ] 1/2 (2.9) 

V T = t ( V TX ) 2 + ( V TY ) 2 ] 1/2 (2.10) 

The line-of-sight angle (a) is defined as 

o= arctan (2.11) 

V -*r ) 

and the rotation rate of the line-of-sight (d) can be 
expressed as 

d = vysin ( y t -o) -l^sintY^-a) (2.12) 

R 

where V T sin(7 X - a) and V M sin(7 M - a) are the target and 
missile velocity components normal to the line-of-sight. 
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III. KINEMATIC REPRESENTATION OP A SKID-TO-TORN MISSILE 

A. INTRODUCTION 

Two basic methods of controlling the attitude of a missile 
to achieve the acceleration commanded by the guidance law are 
skid-to-turn and bank-to-turn. In the skid-to-turn method the 
roll angle is held to a small quantity and is usually 
considered zero in initial calculations. The magnitude and 
orientation of the body acceleration vector is achieved by 
permitting the missile to develop both an angle of attack and 
a sideslip angle. The presence of the sideslip imparts a 
skidding motion to the missile. The bank-to-turn missile 
generally will develop higher lift accelerations than the 
skid-to-turn method, so that the missile requires banking 
maneuvers to properly direct the control vector. To achieve 
the desired orientation, the missile is rolled or banked so 
that the plane of maximum normal force is oriented in the 
desired direction. The present work assumes a skid-to-turn 
missile that is roll stabilized. 

B. ASSUMPTIONS 

A simplified, point mass representation of the missile 
equations of motion will be developed under the following 
assumptions . 
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1. The missile thrust exactly cancels drag. 

2. The orientation of the missile can be described by the 
Euler angles that represent the flight path angles in 
pitch and yaw. 

3 . The seeker head angle rate is a good estimate of the 
line-of-sight rate. 

4. For small angles of attack and sideslip, the velocity 
vector is assumed to be aligned with the body center 
line and if the missile speed is maintained nearly 
constant, the missile acceleration in the x direction 
is fairly close to zero. 

C. SIMPLIFIED EQUATIONS OF MOTION 

Figure 6 illustrates the flight path geometry of the 
missile or target, with the velocity vector aligned with the 
body centerline. From the flight path geometry, the missile 
and target velocity magnitudes are expressed as 

v„ - I J 1/! (3.1) 

for the missile, and 

. t(V„)’.(V„)>.(V„)M lrt (3.2) 

for the target. 

The angles of attack ( 7 _pitch) and sideslip ( 7 _yaw) are 
defined as kinematic quantities depending only on velocity 
ratios. From Figure 6, the missile and target flight path 
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z 




Figure 6. Flight Path Geometry 
angles can be defined as 



Yu pitch * arc tan 



'mz 



/T(vWMvw> 2 >) 



( 3 . 3 ) 



y M yaw * arctan 



( V) 

[ V J 



(3.4) 



Y T pitch * arctan 



rz 






(3.5) 



13 



(3.6) 



Y r yaw * arctan 




If the velocity vector is not aligned with the body 
centerline, the angle of attack is defined as 

Y H pitch » arctan 






(3.7) 



with the sideslip angle in the yaw plane remaining unchanged. 
Figure 7 illustrates the sightline geometry. 




Figure 7. Sight line Geometry 

From the sightline geometry, the line-of-sight angles can be 
determined. 
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a pitch = arc tan 



k JUx t -x m ) 2 *Iy t -y m )*) j 



(3.8) 



a yaw - arctan 




(3.9) 



The range vector (R) can be defined as 



R = /T ( X T -X M ) 2 + ( Y t - Y m ) 2 * ( Z T -Z M ) 2 ) 



(3.10) 



The velocity components can be found from the Euler angle 
transformation. The derivation of the Euler angles can be 
found in Reference 4. The x, y, and z directions are taken as 
the longitudinal, lateral, and normal axis of the missile, 
respectively. The corresponding angles that represent the 
angular displacements are <p (roll) , 9 (pitch) , and i> (yaw) . 
With the missile assumed to be roll stabilized, the roll angle 
will be taken as zero and not included. The Euler 
transformation matrix is given as 



where I points north, J east, and K down. The transformation 
matrix represents the total transformation from the inertial 
coordinate system to the missile body axes. 



'x\ [cos (0) cos (i|r) cos (0) sin(ijr) -sin(0)l[j 
y = -sin(ijf) cos(t(f) 0 J 

z\ |sin(0) cos (i|r) sin(0) sin(i|r) cos(0)J[.k; 
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From the above matrix, the linear components of the 
missiles velocity can be found using the fight path angles in 
the yaw and pitch planes. 

v hx = V„(cos(e)cos(ilr) ) (3.12) 

V = V M (cos ( 0 ) sin (\|r) ) ... 



^ hz ~ "Sih (B) ) (3.14) 

A skid-to-turn missile is controlled by generating the 
required acceleration commands in the pitch and yaw plane. 
From Reference 2, the commanded acceleration in the pitch and 
yaw planes are defined as 

Acom picch = V M pi tch Y M pitch (3.15) 



Acom yav = V M yaw y M yaw (3.16) 

The velocity components of the missile in the pitch and yaw 
planes are defined as 



V H pi tch - V M cos( Y/r Y aw ~ 0 Y aw ) 



(3.17) 



V H yaw = V M cos( pitch ) 



(3.18) 



The autopilot of the missile will convert the commanded 
accelerations into angular rates 7 M pitch and -y M yaw. These 
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angular rates are passed to the control surface servos and are 
converted to the necessary fin deflections required to steer 
the missile to the desired course. 

D. SYSTEM SIGNAL FLOW GRAPH 

Figure 8 shows the total system signal flow graph in the 
y direction, where at launch the missile/target line-of-sight 
is taken as the x axis. The outputs from the seeker is the 
seeker head angle rate, which when multiplied by the 
navigation ratio, becomes the input to the autopilot. The 
autopilot provides the signals required for the missile 
dynamics . 

E. LATERAL AUTOPILOT 

As stated previously, the autopilot receives commands from 
a guidance computer and processes these commands to control 
the deflections of the control surface. There are three 
principle requirements when designing an autopilot. These are, 
quick response, stability, and robustness. The autopilot 
should be able to handle broad variations of the aerodynamic 
parameters. For example, the missile may be required to 
operate over a broad range of Mach numbers and at various 
altitudes which will effect the normal force coefficient, 
which is a function of the Mach number and the altitude. 

In this simulation, a simple lateral autopilot will be 
modeled as a first order lag with a time constant (r) of 0.33. 
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Figure 8. System Flow Graph 
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This has been demonstrated in past simulations to adequate iy 
approximate a more detailed autopilot. 

The lateral autopilot controls the motion of the missile 
in the pitch and yaw planes. As in the case of this study, a 
symmetrical cruciform missile, the pitch and yaw autopilots 
are often identical so only one will need to be modeled. With 
the assumption that the angle of attack is very small, the 
velocity of the missile is aligned with the missile 
centerline. From Figure 8, we have the following expression 

( 3 . 19 ) 

where is the seeker head gimble angle rate and 9 M is the 
angular acceleration of the missile flight path. K is equal to 
( 1/r) . The transfer function of the autopilot is given as 

U (5> . JLM ( 3 . 20 ) 

$ (S) S*K 

F. SEEKER HEAD DEVELOPMENT 

A homing head, when mounted in an airborne missile, is 
called a seeker. The purpose of the seeker is to detect, 
acquire and track a target by sensing radiation or reflected 
energy from the target. In the present study, a narrow field- 
of-view seeker gimballed to the airframe will be designed. The 
seeker maintains the target within this narrow f ield-of-view 
by rotating the platform. 
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If the platform is inertially stabilized, rotation is achieved 
by applying torques which are proportional to the target's 
displacement. From the centerline of the f ield-of-view, the 
equation of motion is 



where T is the applied torque, I is the moment of inertia, and 
/? is the angular acceleration of the seeker head angle. From 
Figure 8 and equation (3.21) the resulting equation of motion 
is 



5 - -f - -K, t P * *i 0 (3,22) 

where K t and K 2 are functions of the time constants of the 
seeker. The transfer function of the seeker head can then be 
expressed as 

P (5) = *i 

o (S) s 2 + kj 5 + K t (3.23) 



G. CONTINUOUS TIME STATE EQUATIONS 
l. Missile Dynamics 

Given the continuous time state equation 
£ ( t) = A x ( t) + B u ( t) 

the missile or target state equation is defined. 
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X 

9 

t 



(3.25) 



X = 



0 1 0 0 0 0 
0 0 0 0 0 0 
0 0 0 1 0 0 
0 0 0 0 0 0 
0 0 0 0 0 1 
0 0 0 0 0 0 





X 




0 0 0 




X 




10 0 




y 




0 0 0 




y 




0 10 




z 




0 0 0 




t. 




■ 

o 

o 

H 



2. Autopilot state Equation 

The pitch and yaw autopilots are identical, so the 
following state equation for the autopilot is defined as 





’-3 0 ' 


Y * Pi t cti 




3 o' 


'0 pitch 


Y = 


. 0 -3. 


Y H yaw 


+ 


0 3. 


0 yaw 



3. Seeker Head State Equation 

The seeker head continuous state equation with a time 
constant of 0.1 seconds in the pitch plane is defined as 



0 1 


'0 pitch 




0 


-100 -20. 


.0 pitch 


+ 


100. 



(3.27) 



The yaw autopilot is identical to the pitch autopilot, so the 
above equation can be easily determined in the yaw plane. 

4. Continuous to Discrata 

Given the continuous time state equations, the 
discrete time equations are defined. 
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x (k + 1) = 4> x(k) + T u(k) 



(3.28) 



y (k + 1) = C x(k) + I? u(k) (3.29) 

The simulation study will use the discrete time representation 
of the given continuous time equations. 

H. SIMULATION 

This section presents the results of the missile/target 
engagement using the proportional navigation scheme. The 
following assumptions are made: 

1. The missile is limited to 20 g's in either the yaw or 
pitch plane. 

2. The target is capable of a 5 g maneuver. 

3. The seeker head system is noise free. This will form a 
basis from which the following chapters can be referred 
to. 

4. The minimum difference between the targets position and 
the missiles position, will be used as an estimate of 
the miss distance parameter. 

5. The origin is taken as coordinates (0,0,0) in the x,y,z 
plane. 

l. Constant Velocity Target 

The first scenario will be a constant velocity target 
with the following initial conditions. 
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v mx(°) = 2000 feet/sec 
X T (0) = 20,000 feet 
Z T (0) = 50 feet 
v tx(0) = 1000 feet/sec 

All other initial conditions are zero. As shown in Figure 9, 
a successful intercept occurred for the constant velocity 
target. The minimum distance between the target and missile 
was 0.7665 feet. Figure 10 is a plot of the missile/target 
trajectories looking down on the X-Y plane. It shows the 
target starting at 20,000 feet on the X axis with the missile 
starting at the origin. Figure 11 is a plot of the 
missile/target trajectories in the Y-Z plane with the missile 
starting at the origin. Figure 12 is a plot of the missile 
commanded acceleration parameter. The acceleration increases 
rapidly during the initial phase of the engagement as force is 
applied to the missile and then drops off until just prior to 
intercept. Figure 13 shows the seeker head angle rate in the 
yaw plane and Figure 14 is a plot of the seeker head angle 
rate in the pitch plane. These parameters also increase 
rapidly during the first few seconds as the target is starting 
out on its flight profile and then decreases in time as the 
missile is on course due to the acceleration commands. 
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Figure 9. Range vs Time 



MISSILE-TARGET TRAJECTORY IN THE X Y PLANE 




Figure 10. Trajectory in the X-Y Plane 
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MISSILE. TARGET TRAJECTORY IN THE YZ PLANE 




Figure 11. Trajectory in the Y-Z Plane 
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Figure 12. Commanded Acceleration 
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Figure 13. Seeker Head Angle Rate in the Yaw Plane 
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Figure 14. Seeker Head Angle Rate in the Pitch Plane 
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2. Constant Acceleration Target 

The second scenario will be a constant acceleration 
target with the following initial conditions: 

Vmx(°) = 2000 feet/sec 

x T (0) = 20,000 feet 

Z T (0) = 50 feet 

V TY (0) = 500 feet/sec 

a tx(°) = -5.0*32.2 f eet/sec~2 

a ty(°) = 3.5*32.2 feet/sec A 2 

a tz(°) = 3.5*32.2 feet/sec*2 

All other initial conditions are zero. Figure 15 is a plot of 
the minimum range between the missile and target. The minimum 
distance was 1.13 feet. This was as expected due to the 
changing target trajectory. Figure 16 is the plot of the 
missile/target trajectory in the X-Y plane and Figure 17 is 
the missile/target trajectory in the X-Z plane. Both plots 
demonstrate the ability of the missile using the proportional 
navigation guidance law, to effectively intercept the target. 
Figure 18 is a plot of the missile acceleration vs time. As 
expected, in contrast to the constant velocity simulation, the 
initial portion of the engagement is characterized by a low 
acceleration magnitude. As the engagement progress, however, 
due to the missile having to pull more g's in order to 
compensate for the more rapid trajectory of the target, the 
missile acceleration increased rapidly. 
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Figure 19 is a plot of the targets acceleration vs time. 
Figures 20 and 21 are the seeker head angle rates in the yaw 
and pitch planes, respectively. As with the acceleration of 
the missile, similar characteristics were realized for these 
parameters due to the fact that the missile is constantly 
changing its flight path to compensate for the acceleration of 
the target. 
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Figure 15. Range vs Time 



28 





30 



MISSILE-TARGET TRAJECTORY IN THE X-Y PLANE 




Figure 16. Trajectory in the X-Y Plane 



MISSILE-TARGET TRAJECTORY IN THE X-Z PLANE 




Figure 17. Trajectory in the X-Z Plane 
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Figure 20. Seeker Head Angle Rate in the Yaw Plane 
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Figure 21. Seeker Head Angle Rate in the Pitch Plane 
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IV. LUENBERGER OBSERVER 



A. INTRODUCTION 

The purpose of a Luenberger observer for the proportional 
navigation model developed in Chapter III, is to estimate the 
variables in the missile/target engagement that might not be 
measured directly or to reduce noise on the states that are 
measured by the missile. It is not practical to design a 
missile to determine all inputs. Such inputs could be evasive 
maneuvers of the target, initial conditions, noise level, or 
the jamming techniques employed. The use of a Luenberger 
observer is beneficial in these situations. 

B. OBSERVER DESIGN 

The seeker head angle rate is assumed to be a good 
approximation of the line-of-sight rate. The observer will be 
designed to estimate the states of the seeker head, since the 
actual line-of-sight rate of the seeker head can be 
contaminated by noise or jamming. The following equations are 
described in Reference 4. Given a system 

x(k + 1) = <bx(k) + Au(Jc) (4.1) 

and the output equation 

y(k) = C x(k) + D u(k) (4.2) 

The observer will approximate x(k) with inputs y(k) and u(k). 
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The observer equation is given as 

X(k * 1) = F k(k) + G y(k) + H u(k) (4.3) 

with xhat(k) being the estimate of the original system. The 
error between the state x(k) and the estimated state can be 
defined as 



e(k * 1) = x(k + 1) - k{k * l) (4.4) 



e(k * 1) * #x(k) +A u(k) - F k{k) - G y(k) - H u(k) 



(4.5) 



which results in the error equation of the form 

e{k * 1) = (4>-GC)x(Jc) - F k(k) + (A -GD-H)u{k) (4.6) 

If we select F as 



F = <b-GC 



(4.7) 



Equation (4.4) becomes a state equation of the following form 



e(Jc + 1) = {<b-GC) e(k) + ( A -GD-H) u ( k) 



(4.8) 



By selecting H as 



H = A -GD (4.9) 

Equation (4.6), with input u(k), results in a state equation 
without input, of the following form 



e(k + 1) = (4>-GC) e(k) (4.10) 

If the original system is observable, by selecting the 
observer gain G, we can design the observer such that the 
error approaches zero at a rate that is desired based on the 
design characteristics. 
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The unforced system 



x(k ♦ 1) = $ x (4.11) 

with the observation equation 

y(k) = C x(k) (4.12) 

is observable if and only if the rank (N) of the observability 
test matrix 



N = [C' <& / C / (<X>) Jc - l C / ] (4.13) 

is equal to k, the order of the system [Ref. 5]. 

The observer design begins by converting the state 
equation for the seeker head, Equation (3.27) , from continuous 
time to discrete time representation. With a sampling time of 
0.05 and an unknown disturbance w(k) the equation is given as 



* (Jc*l) 
P(JC-I) 



0.9098 


0.0303 


'P(Jc) 




0.0902' 


a(k) ♦ 


0.0013 


-3.0327 


0.3033, 


P(lc). 




3.0327, 


.0.0513, 



(4.14) 



with the output equation measuring the seeker head angle rate 
given as 



y(k) 



[o 1] 



P(k) 

0(ic) 



(4.15) 



The system is observable, so the observer gain can be 
determined. It is usually desirable to select the eigenvalues 
of the observer so that the error between the state and 
estimated state will rapidly become small and will not be 
highly responsive to noise. Because the system is observable, 
the eigenvalues of the observer can be arbitrarily positioned. 
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If the observer gain is high, the observer will quickly drive 
the initial condition error to zero but will make the observer 
more susceptible to noise. By choosing a low observer gain, 
the observer response will be slower but the observer will be 
less affected by noise or jamming. Previous studies have shown 
that a good compromise between the two would be to choose 
continuous time poles of P[ and p 2 to be between the values of 
8.0 to 12.0. Continuous time poles of 10.5 were chosen for the 
observer which results in the discrete time poles given as 

Poi = P02 = e ( ' 10,5 * °-° 5) =0.5916 (4.16) 
The observer gain matrix is then given as 



-0.0031 

0.0299 



(4.17) 



The F matrix is equal to (* - GC) and is given as 



‘0.9098 


0.0303' 




-0 .0031 


-3 .0327 


0.3033. 




. 0.0299 . 



[0 1 ] 



(4.18) 



The H matrix is equal to (A - GD) and with D equal to zero the 
H matrix reduces to a. The observer equation for the seeker 
head is given as 



[?<*♦!>) 



0.9098 0 . 0297l.£ 


-0 . 003ll , f% 


0.0902 


[-3 .0327 0 . 3094f ^ 


[o.0299 F (JC) * 


3.0327 



oik) 



(4.19) 
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C. SIMULATION RESULTS 

The two previous scenarios given in Chapter III were used 
to test the observer design. Random noise was injected into 
the seeker head to determine if the Luenberger Observer would 
reduce the noise of the states that were directly measured by 
the seeker head and generate the required signals necessary 
for the autopilot. 

l. Constant Velocity Target 

The same initial conditions presented in Chapter III 
for the constant velocity engagement will be used for 
comparison in this section. Figures 22 through 27 are plots of 
the constant velocity engagement without the Luenberger 
Observer but with a random disturbance. Figures 22 and 23 are 
the plots of the trajectories in the X-Y and X-Z planes, 
respectively. The miss distance was 15.6 feet. The trajectory 
in the Y-Z plane depicts how the noise input to the seeker 
head adversely effects the missiles trajectory. Figures 24 and 
25 are the seeker head angle rates in the yaw and pitch plane, 
respectively and Figure 26 is the plot of the commanded 
acceleration as a result of the seeker head angle rates being 
subjected to noise. Figure 27 is the plot of the random noise 
generated by the program. 
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MISSILE- TARGET TRAJECTORY IN THE X-Y PLANE 
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Figure 22. 


Trajectory in the X-Y Plane With Noise 
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MISSILE- TARGET TRAJECTORY IN THE Y-Z PLANE 
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Figure 23. Trajectory in the Y-Z Plane with Noise 
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Figure 24. Seeker Head Angle Rate in the Yaw Plane 
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Figure 25. Seeker Head Angle Rate in the Pitch Plane 
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Figure 26. Conunanded Acceleration with Noise 




Figure 27. Random Noise Generated vs Time 
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The following plots depict the results with the 
observer included in the simulation with the same initial 
conditions. Figures 28 and 29 are the trajectories in the X-Y 
and Y-Z planes, respectively. In contrast with the Y-Z plane 
plot without the observer, the missiles trajectory is less 
effected by the random disturbance. A miss distance of 1.56 
feet occurred. Figures 30 and 31 are the plots of the seeker 
head angle rates in the yaw and pitch plane, respectively. 
When compared to the angle rate plots of Chapter III, it is 
clear that the observer has reduced the effect of noise to the 
point that the plots are very similar. Figures 32 and 33 are 
the plots of the commanded acceleration and the random noise 
generated, respectfully. The commanded acceleration is also 
less effected by the random disturbance when the observer is 
made part of the missile. 
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Figure 28. Trajectory in the X-Y Plane with Observer 




Figure 29. Trajectory in the Y-Z Plane with Observer 
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Figure 30. Seeker Head Angle Rate in the Yaw Plane 
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Figure 31. Seeker Head Angle Rate in the Pitch Plane 
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Figure 32. Commanded Acceleration with Observer 




Figure 33. Random Noise Generated vs Time 
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2. Constant Acceleration Target 

The results for the constant acceleration engagement 
were as expected, with the observer reducing the effects of 
noise so that intercept occurred. Figures 34 through 39 are 
plots of the constant acceleration model without the observer. 
Figures 34 and 35 are plots of the trajectories in the X-Y and 
X-Z planes, respectively. These plots do not show the 
variations in the trajectories due to the magnitudes of the 
scales. Figures 36 and 37 are plots of the seeker head angle 
rates in the yaw and pitch planes, respectively. When compared 
to the plots in Chapter III, it is obvious that the missiles 
performance is degraded. A miss distance of 26.5 feet 
occurred. Figures 38 and 39 are plots of the commanded 
acceleration and the random noise generated. 
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MISSILE-TARGET TRAJECTORY IN THE X-Y PLANE 




Figure 34. Trajectory in the X-Y Plane with Noise 



MISSILE- TARGET TRAJECTORY IN THE X-Z PLANE 




Figure 35. Trajectory in the X-Z Plane with Noise 
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Figure 36. Seeker Head Angle Rate in the Yaw Plane 
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Figure 37. Seeker Head Angle Rate in the Pitch Plane 
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Figure 38. Commanded Acceleration with Noise 
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Figure 39. Random Noise Generated vs Time 
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The following plots depict the constant acceleration 
engagement with the Luenberger Observer. The observer reduced 
the noise so that intercept occurred with a miss distance of 
4.5 feet. The plots in the X-Y and X-Z planes are not included 
due to the magnitude of the scales. Figures 40 and 41 are the 
plots of the seeker head angle rates in the yaw and pitch 
planes, respectfully. When compared with the constant 
acceleration model in Chapter III it is obvious that the 
observer is beneficial in reducing the random disturbance 
added to the seeker head. Figure 42 is the commanded 
acceleration of the missile. 

In the present day environment, with fast maneuvering 
targets, the simulations in this chapter demonstrate the need 
for the observer concept to be included in the missile design 
and also in computer simulation studies. 
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Figure 40. Seeker Head Angle Rate in the Yaw Plane 
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Figure 41. Seeker head Angle Rate in the Pitch Plane 
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Figure 42. Commanded Acceleration with Observer 

The computer program used to generate the plots in this 
chapter is included in the appendix, of which a portion is the 
work completed previously by LT F. Lukenbill in our missile 
and avionics class. 
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V. GROUND TARGET OBSERVER 



A. INTRODUCTION 

Air defense is becoming significantly more complex due to 
the increased speed and manueverability of hostile aircraft 
and the enemy's use of electronic countermeasures (ECM) . 
Attacking aircraft will employ jamming as well as other forms 
of ECM to degrade the tracking performance of surface-to-air 
missiles. Therefore, many modern surface-to-air missile 
systems have target observers that track incoming aircraft 
from the ground and uplink the deviations in target position 
and velocity to the missile. 

In this chapter, a simplified Kalman Tracking Filter is 
developed. The tracker consists of velocity and position 
estimation. Once the missile is launched, target deviations 
are computed by subtracting the target observer estimates from 
the current radar estimates. At specified intervals in time, 
the errors in position and velocity are uplinked to the 
missile. After reception, the missile converts the errors into 
the missile body reference frame via the attitude Euler angles 
and then into the line-of-sight frame using the seeker head 
angle transformation. The estimated angular rate deviations 
are then added to the missile's current estimate and the 
appropriate acceleration commands are generated. 
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B. KALMAN FILTER 

The purpose of the Kalman Filter is to keep track of the 
states of a system given a sequence of noisy measurements. 
This is accomplished by recursively updating an estimate of 
the state by processing a sequence of noisy observations so 
that the effect of the measurement errors are reduced as much 
as possible. 

First, previous estimates of the state x, and the 
covariance matrix P, are extrapolated one step ahead based on 
the assumed systems dynamics. These values are used to compute 
a set of optimum weights or Kalman Gains. The gains are 
applied to the prediction and to a new observation which 
provides an updated estimate of the state and the covariance 
matrix. 

It is assumed that the target tracker has no a priori 
information on the target's position or velocity. Also, the 
error covariance matrix will be set sufficiently large in 
order that the filter gain will be able to make the estimate 
converge to the true value and make the simulation independent 
of any initial estimate states. The random forcing function W k 
is zero mean with covariance Q k and the measurement noise V k is 
zero mean with covariance Rfc. The Kalman Filter algorithm is 
given in Equations (5.1) through (5.5) in the order that the 
filter is implemented in the computer program. The algorithm 
is given as 
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g* - *Wi, * r w Pi * r ♦ *) -i 



(5.1) 




(5.2) 



^(Jc|Jc) “ (I ~ H) Pjjfijf.D 



(5.3) 




(5.4) 



^(Jc*l|Jc) = ® ^(Jc|Jc) + 0 



(5.5) 



From Reference 2, the following covariance matrices are 
defined. The purpose of the random forcing function covariance 
matrix Q, is to account for model inaccuracies or for a target 
that has maneuvered. The covariance matrix Q is given as 



Since the actual range of the target is unknown, an estimate 
of the range is used in determining R*. The covariance of 
measurement error matrix will then be based on a three degree 
field of view radar, observing the target approaching an area 
20,000 feet from the ground based tracker. 



0 0 fi60 2 0 0 



Q » 0 0 » 0 160 2 0 

0 0 o\ L 0 0 160 2 



(5.6) 
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500 2 0 0 

0 125 2 0 (5.7) 

0 0 12 . 5 2 



The measurement equation y k is given as 

y(k) = H x(k) + v ( k) (5.8) 

where H is the measurement matrix, a linear function of the 
state vector, and u(k) the measurement noise. 

C. UPLINK 

The uplink from the ground target observer consists of 
deviations in target position and velocity. The deviations are 
calculated by determining the difference between the actual 
noisy measurement and the predicted state estimate. 



y(k) - H •£(*(*_!) 



(5.9) 



Upon reception, the missile transforms the deviations into the 
missile body frame. 







^by 


II 


Ybz. 


. 



cos (0) cos (ijr) cos (0) sin(\|r) -sin(0) 
-sin(i|i) cos(\|r) 0 

sin(0) cos (i|r) sin(0) sin(ijr) cos(0) 





Vx 




Vy 


■ 


y z . 



(5.10) 



where V x , V Y , and V z are the velocity deviations in the x, y, 
and z directions, respecttively , and Vbx, V by , and V BZ are the 
transformed velocities in the missile reference frame. The 
range deviations are calculated the same and are defined! 
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(5.11) 



■KflX 




cos (0) cos (l|/) 


cos (0) sin (\|r) -sin(0) 




R x 


r by 


= 


-sin(ijf) 


cos(ijr) 0 




R r 


r bz. 




sin (0) cos (i|f ) 


sin (0) sin (i(r) cos(0) 




Rz. 



where R x , R Y/ and R z are the range deviations and R BX , R BY and 
R bz are the transformed range deviations in the x, y, and z 
directions, respectively. The transformation to the line-of- 
sight reference frame is accomplished by using the seeker head 
angle transformation and is given as 







cos (o p ) cos (o r ) 


cos (<j p ) sin(o r ) -sin(Op) 




' v ax 


Ksy 
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-sin(a y ) 


cos(o r ) 0 




V.r 


V 3Z. 




sintap) cos (o r ) 


sin(Op) sin(o Y ) cos(o p ) 







where a ? and cr Y are the line-of-sight angles in the pitch and 
yaw plane, respectively. V^, V SY , and V sz are the transformed 
velocities. The range is transformed using the same 
transformation matrix and is given as 



R sx 




COS ( Op) cos (o y ) 


cos (o p ) 8in(o r ) -sin(o p ) 




Rax 


Rsy 


M 


-sin(a y ) 


cos(a r ) 0 




Rby 






sin (Op) cos (o r ) 


sin(Op) sin(o r ) cos (o p ) 




Rbz. 



where R^, Rg Y , and Rj Z are the transformed ranges in the x, y, 
and z directions, respectively. The final process in the 
uplink transformation is to generate the angular rate 
deviations. These are calculated by dividing the velocity 
deviations by the range and the magnitude of the range 
deviations. The line-of-sight rate in the yaw plane is given 
as Equation (5.14). 
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(5.14) 



d 




where R D is given as 



r d = yr r 2 sx * rIy ♦ r 2 sz ) 



(5.15) 



and the line-of-sight rate in the pitch plane defined as 



D. SIMULATION OF THE TARGET OBSERVER 

The constant velocity target in a noisy environment will 
be the first test case for the ground target observer. Figures 
43 through 51 summarize the extensive estimator performance 
calculations for the constant velocity scenario. Figures 43 
through 48 are the histories of the six elements of the state 
vector. An inspection of these figures indicates that the 
estimator tracks the system in position and velocity quite 
well. Figures 49 through 51 are examplese of the Kalman Filter 
Gains vs time. Figures 52 through 54 are the response of the 
missile to the uplink commands. The deviations in position and 
velocity were uplinked at one second intervals. The miss 
distance parameter was observed to be 0.44 feet. 




(5.16) 
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Figure 43. 


Target Position in the X Direction 
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Figure 44. Target Velocity in the X Direction 
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Figure 45. 


Target Position in the Y Direction 
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Figure 46. Target Velocity in the Y Direction 
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Target Position in the Z Direction 
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Figure 48. Target Velocity in the Z Direction 
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Figure 49. 


Kalman Gains in the X Direction 
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Figure 50. Kalman Gains in the Y Direction 



60 







KALMAN GAINS u the Z DIRECTION 








20 














18 










- 




16 


- 












14 




(POSITION GAIN) 










12 




t + + (VELOCITY GAIN) 






- 




10 














8 


- 




* 








6 


- 




• 








4 






. 










. 












2 


_ ♦ 






r 








"... V , 


























°« 


1 


2 3 4 5 


6 


7 


I 


1 






TIME 









Figure 51. Kalman Gains in the Z Direction 
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Figure 52. Seeker Head Angle Rate in the Yaw Plane 
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Figure 53. Seeker Head Angle Rate in the Pitch Plane 
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Figure 54. Acceleration Commanded 
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For the case of a maneuvering target, the constant 
acceleration model was used to determine filter performance. 
As expected, filter performance was marginal. Filter 
convergence to an accuracy attained in the constant velocity 
target did not occur. Figures 55 through 63 summarize the 
estimator performance for the constant acceleration scenario. 
Figures 55 through 60 are the six elements of the state 
vector. As shown in these plots, the estimator does not track 
the maneuvering target as well as the constant velocity model. 
Figures 61 through 63 are the Kalman Filter Gains generated 
and Figures 64 through 66 are the missile's response to the 
uplink commands. The uplinks were transmitted using the same 
time interval as in the constant velocity scenario. A miss 
distance of 8.6 feet occurred. 
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Figure 55. 
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Figure 56, 


, Target Velocity in the X 


Direction 
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Figure 57. Target Position in the Y Direction 
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Figure 58. Target Velocity in the Y Direction 
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Figure 59. Target Position in the Z Direction 
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Figure 60. Target Velocity in the Z Direction 
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Figure 61. Kalman Gains in the X Direction 
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Figure 62. Kalman Gains in the Y Direction 
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Figure 63. Kalman Gains in the Z Direction 
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Figure 64. Seeker Head Angle Rate in the Yaw Plane 
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Figure 65. Seeker Head Angle Rate in the Pitch Plane 
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Figure 66. Commanded Acceleration 
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The computer program used to generate the plots in this 
chapter is included in the appendix. 
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VI. CONCLUSIONS AND RECOMMENDATIONS 



A. CONCLUSIONS 

The use of a three-dimensional model in simulation studies 
is beneficial in that the trajectories of the missile and 
target can be displayed in all three axes, therefore, a better 
understanding of the behavior of the systems involved. 

The simulation studies in Chapter III demonstrated that 
the Proportional Navigation guidance law effectively achieves 
intercept of the target for both the maneuvering and non- 
manuvering target. Chapter IV demonstrated the need for the 
observer concept to be utilized in missile design to reduce 
the effects of noise. Chapter V showed an additional means of 
countering jamming or ECM by employing the ground target 
observer concept. This demonstrated that the uplinks in target 
position and velocity, for the constant velocity model, were 
beneficial in a noisy environment. 

B. RECOMMENDATIONS 

Further studies could include a Kalman Filter in the 
missile design instead of the Luenberger Observer. It is 
recommended that future studies include the use of a two part 
Kalman Filter with the ground target observer. One that 
estimates the states and the other being a maneuver detector. 
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The results could be compared with the present filter along 
with ECM effects to determine the effectiveness of the two. 
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APPENDIX 



A. MISSILE-TARGET ENGAGEMENT PROGRAM 



% 

% This program simulates a 3 dimensional target/missile 

% engagement using basic proportional navigation with 

% an asymptotic observer. 



% DEFINE STATES 

% 

% Missile 

% 



% 


ms=[xm 


- missile 


% 


xdm 


missile 


% 


ym 


- missile 


% 


ydm 


- missile 


% 


zm 


- missile 


% 


zdm 


- missile 



% 

Am=[ 0 10000 
0 0 0 0 0 0 
0 0 0 1 0 0 
0 0 0 0 0 0 
0 0 0 0 0 1 
000000 ]; 



position in x direction 
velocity in x direction 
position in y direction 
velocity in y direction 
position in z direction 
velocity in z direction] 



Bm=[0 0 0 
10 0 
0 0 0 
0 10 
0 0 0 

0 0 1]; 

% 

% Seeker Head 

% 

% betap=[beta pitch - seeker head pitch angle 

% betad pitch - seeker head pitch angle rate 

% betay« [beta yaw - seeker head yaw angle 

% betad yaw - seeker head yaw angle rate] 

% 

% 



Asp= [ 0 1;-100 -20.0]; 
Asy= [ 0 1 ; -100 -20.0]; 
Bsp=*[0;100] ; 
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Bsy=[0;100] ; 

% 

% Autopilot 

% 

% gammamp= [ gammam pitch - pitch body angle 

% gammadm pitch - pitch body angle rate ] 

% gammamy= [ gammam yaw - yaw body angle 

% gammadm yaw - yaw body angle rate] 

% 

Ap=[-3 0 ; 0 -3]; 

Bp= [3 0 ; 0 3 ] ; 

% 

% Target 

% 



% 


ts=[xt 


target 


position 


in 


X 


direction 


% 


xdt 


target 


velocity 


in 


X 


direction 


% 


yt 


target 


position 


in 


y 


direction 


% 


ydt 


target 


velocity 


in 


y 


direction 


% 


zt 


target 


position 


in 


z 


direction 


% 


zdt 


target 


velocity 


in 


z 


direction] 



% 

At= [ 0 10000 
0 0 0 0 0 0 
0 0 0 1 0 0 
0 0 0 0 0 0 
0 0 0 0 0 1 
000000 ]; 

Bt= [ 0 0 0 
10 0 
0 0 0 
0 10 
0 0 0 
ooi]; 

% 

% DISCRETE REPRESENTATION 

% 

dt= . 05 ; 

[phisp,delsp]= c2d(Asp,Bsp,dt) ; 
[phisy,delsy]= c2d(Asy, Bsy ,dt) ; 
[phim,delm]= c2d(Am,Bm,dt) ; 
[phia,dela]= c2d (Ap, Bp, dt) ; 
[phit,deltj= c2d(At,Bt,dt) ; 

tfinal= 20.0; 
kmax=tf inal/dt +1; 

% INITIAL VALUES 



74 



% 

% navigation ratio 

% 

NR= [4.0 0 ; 0 4.0]; 



% initial seeker head angles and angle rates 

% 

beta_pitchO =0.0; 
betad_pitch0=0 . 0 ; 
beta_yaw0 =0.0; 
betad_yawO =0.0; 

betap ( : , 1 ) = [ beta_pitchO 

betad_pitchO ] ; 
betay ( : , 1 ) = [ beta_yaw0 

betad_yawO ] ; 

% initial observer angles and angle rates 

betae_pitch0=0 . 0 ; 
betade_pitch0=0 . 0 ; 
betae_yaw0=0 . 0 ; 
betaed_yaw0=0 . 0 ; 

% 

% 

betaep( : , 1) =[betae_pitchO 

betade_pitchO] ; 
be taey ( : , 1 ) = [ betae_yawO 

betaed_yawO] ; 

% 

% 

% initial missile body angle rates 

% 

gammadm_pitch0=0 . 0 ; 
gammadm_yawO =0.0; 

gammadm (:,!)=[ gammadm_pitchO 

gammadm_yawO ] ; 

% 

% initial states for missile flight profile 

% 

xmO =000.0; 
xdm0=2000.00; 
ymO =000.0; 
ydm0=00; 
zmO =00; 
zdm0=00 ; 
ms ( : , 1) =[xm0 
xdmO 
ymO 
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ydmO 

zmO 

zdmO ] ; 

% 

% initial states for target flight profile 

% 

xto =20000; 

xdt0=0000 ; 
ytO =050.0; 
ydt0=0300 . 0 ; 
z to =0051.0; 
zdt0=000 . 0 ; 

ts( : , 1) =[xtO 
xdtO 
ytO 
ydtO 
zto 

zdtO] ; 

% 

% initial range information 

% 

rxO= (xtO-xmO) ; 
rx(l)=rxO; 
ryO= (ytO-ymO) ; 
ry (l) =ryO; 
rzO= (ztO-zmO) ; 
rz ( 1) =rzO; 

rmO=sqrt (xmO*2 + ymO*2 + zm0*2) ; 
rm ( 1 ) =rmO ; 

rtO=sqrt (xtO A 2 + yt0~2 + zt0*2) ; 
rt(l)=rtO; 

rO =sqrt ( (xtO-xmO) ~2 + (ytO-ymO) *2 + (ztO-zmO) ~2) ; 

r ( 1 ) =r 0 ; 



dist=[ 0.0013; 0.0513] ; 
rand( 'normal' ) ; 

% 

% initial time 

% 

time (l) =0.0; 

% 

% SIMULATION 

% 

for (i=l:kmax-l) 

% missile and target flight path angles 

gammam_pitch( i) =atan2 (ms (6, i) , sqrt (ms (2, i) *2+ms (4 , i) *2) ) ; 
gammam_yaw(i)=atan2(ms(4,i) ,ms(2, i) ) ; 
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gammat_pitch(i)=atan2 (ts(6,i) , sqrt (ts (2 , i) A 2+ts(4,i) A 2) ) ; 
gammat_yaw(i) =atan2 (ts (4 , i) , ts(2 , i) ) ; 



% target velocity (magnitude) 

vt ( i) =sqrt (ts (2 , i) A 2+ts (4 , i) A 2+ts (6 , i) A 2 ) ; 
vm (i) =sqrt (ms (2 , i) A 2+ms (4 , i) A 2+ms (6 , i) A 2 ) ; 

% 

% line-of-sight angles 

% 

sigma_pitch(i) =atan2 ( (ts(5, i) -ms (5, i) ) , sqrt ( (ts(l, i) -ms(l, i) 
) * 2 • 

+ ( ts ( 3 , i ) -ms ( 3 , i ) ) A 2 ) ) ; 



sigma_yaw(i) =atan2 ( (ts (3 , i) -ms (3 , i) ) , ( (ts (1, i) -ms (1 , i) ) ) ) ; 

sigma ( : , i) =[ sigma_pitch (i) 

sigma_yaw(i) ]; 



% update seeker head 

% 

wk(i) =rand* . 9 ; 



betap( : , i + l)=phisp*betap( : , i) + 
delsp*sigma_pitch (i) +dist*wk ( i) ; 

betay ( : , i+1) =phisy*betay ( : , i) + delsy*sigma_yaw(i) +dist*wk(i) ; 
% 

% set up observer estimates 

% 

yO(l) =0.0; 
yl(l)=0.0; 

% 

% 

gain* (-0.0031; 0.0299] ; 
c-(0 1]; 
f*phisp-gain*c; 
yO ( i+1) *c*betap ( : , i) ; 

betaep( : , i+1) =f *betaep ( : , i) + gain*y0(i) + 

delsp*sigma_pitch(i) ; 

yl (i+1) =c*betay ( : , i) ; 
f l=phisy-gain*c; 

betaey ( : , i+1) =f l*betaey ( : , i) + gain*yl(i) + 

delsy*sigma_yaw(i) ; 
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% 

% 

% seeker head angle rate vector 

% 

betad ( : , i ) = [ betaep ( 2 , i ) 

betaey (2 , i) ] ; 

% 

% update autopilot 

% acceleration commanded = vm(i) * gammadm(i) 

% gammadm(i) = Nav constant * line of sight rate 

gammadm( : , i+1) =phia*gammadm ( : , i) + dela*NR*betad(:,i); 

vm_pitch(i)=vm(i) *cos(gammam yaw(i) -sigma_yaw(i) ) ; 
vm_y aw ( i ) =vm ( i ) *cos ( gammam_pTtch ( i ) ) ; 



% limit commanded acceleration to approximately 20 g's 

% 

if (vm_yaw(i) *gammadm(2 , i) ) <=640 . 0 

amc_y (i) = (vm_yaw(i) *gammadm(2 , i) ) ; 

else 

amc_y(i) =640. 0; 

end 

% 

% 

if (vm_pitch(i) *gammadm(l, i) ) <=640. 0 

amc_p(i) =(vm_pitch(i) *gammadm(l, i) ) ; 

else 

amc_p(i) =640. 0; 

end 

% 

% magnitude of commanded acceleration 

acom(i)=sqrt (amc_y(i) ~2 + amc_p(i)*2); 

% missile acceleration vector components 



xddm(i)=(gammadm(l / i) *ms(6, i) *cos ( gammam_y aw ( i ) ) ) ; 
yddm(i) =amc_y (i) *cos (gammam_yaw(i) ) ; 
zddm(i)=amc_p(i) *cos(gammam_pitch(i) ) ; 

% total missile acceleration magnitude 

% 

am(i)=sqrt (xddm(i) A 2 + zddm(i)~2 + yddm(i) A 2 ); 

% 

% missile input vector 
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<*><*><*><*> <#><#><#> <#><#><#> <#><#><#> <#><#> <#> <#><#><#> 



% 



um= [ xddm ( i ) 
yddm ( i ) 
zddm(i) ] ; 

update missile 

ms ( : , i+1) =phim*ms ( : , i) + delm*um; 
set target acceleration components 
if (r (i) <=21000 . 0) 

xddt(i)=-5.0*32.2*sin(sigma_yaw(i) ) ; 
yddt(i)* 3.5*32.2*cos(sigma_yaw(i) ) ; 
zddt(i)* 3.5*32.2*cos(sigma_pitch(i) ) ; 

else 

xddt(i)=0.0; 
yddt (i) =0 . 0; 
zddt ( i) =0 . 0 ; 

end 



target acceleration magnitude 

at(i)=sqrt(xddt(i) A 2 + yddt(i) A 2 + zddt(i) A 2); 

target input vector 

ut=[xddt (i) 
yddt (i) 
zddt (i) ] ; 

update target 

ts ( : , i+1) »phit*ts ( : , i) + delt*ut; 



missile and target trajectory data 

missile(i, : ) »[ms (1, i) ms(3,i) ms(5,i)]; 
target(i,:) =(ts(l,i) ts(3,i) ts(5,i)]# 

% 

% update range information 
% 



r ( i+1) «sqrt ( (ts (1, i+1) -ms (1, i+1) ) A 2 + 

(ts (3, i+1) -ms (3, i+1) ) *2. . . 

+ (ts(5, i+1) -ms (5, i+1) ) A 2) ; 
rm(i+l) =sqrt (ms (1, i+1) A 2 + ms (3 , i+1) A 2 + ms (5, i+1) A 2) ; 
rt(i+l)=sqrt(ts(l,i+l) A 2 + ts(3,i+l) A 2 + ts(5,i+l) A 2) ; 
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% update time vector 

% 

time ( i+1) =time ( i) + dt; 

% 

% check to see if engagement is at closest point of 

approach 

% 

if(r(i)<r(i+l)) , break, end 

end; 

% 

% PLOT RESULTS 

% 

% range data 

% 

plot (time, r) , grid, xlabel ( 'TIME' ) , y label ( 'FEET' ) 
title ('RANGE VS TIME') 

text( 7,20000, [ 'rng=' ,num2str(r(i) ), ' feet']); 
meta msltgt, pause, clg 

subplot (211) ,plot (time,rm) , grid, xlabel ( 'TIME' ) ,ylabel( ' FEET' ) 
title ( 'MISSILE RANGE VS TIME'); 

subplot (212 ) , plot (time, rt) , grid, xlabel ( 'TIME' ) ,ylabel( 'FEET' ) 
title ( 'TARGET RANGE VS TIME'); 
meta jerryl, pause, clg 
% 

% 

plot (ms ( 1 , : ) , ms ( 3 , : ) , ts ( 1 , : ) , ts ( 3 , : ) ) , grid 

title ( 'MISSILE-TARGET TRAJECTORY IN THE X-Y PLANE') 

xlabel ('X AXIS IN FEET' ) , ylabel ( ' Y AXIS IN FEET') 

gtext ( 'MISSILE' ) 

gtext( 'TARGET' ) 

meta msltgt, pause, clg 

plot (ms ( 1, : ) ,ms (5, : ) , ts (1, : ) , ts (5, : ) ) ,grid 

title ( 'MISSILE-TARGET TRAJECTORY IN THE X-Z PLANE') 

xlabel ('X AXIS IN FEET' ), ylabel (' Z AXIS IN FEET') 

gtext ('MISSILE') 

gtext ( 'TARGET' ) 

meta msltgt, pause, clg 

% missile and target velocity 

% 

plot (time ( 1 : i) ,vm) , grid, xlabel ( 'TIME' ) , ylabel (' FEET/ SEC' ) 
title ( 'MISSILE VELOCITY VS TIME'); 
meta jerryl, pause, clg 

plot (time ( 1 : i) ,vt) , grid, xlabel ( 'TIME' ) , ylabel (' FEET/ SEC' ) 
title( 'TARGET VELOCITY VS TIME'); 
meta jerryl, pause, clg 

% missile and target acceleration 
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plot (time (1 : i) , am) , grid, xlabel ( 'time' ) , y label ( ' feet / sec A 2 ' ) 
title( 'Missile acceleration vs Time') 
meta jerryl, pause ,clg 

plot (time ( 1 : i) , at) , grid, xlabel ( ' time ' ) , y label ( ' feet / sec A 2 ' ) 
title ( 'Target acceleration vs Time') 
meta jerryl , pause , clg 

% seeker head angle rates 

plot (time ( 1: i-1) , betay (2 , (1: i-1) ) ) , grid, xlabel ( ' time ' ) ,ylabe 
1 ( ' radians ' ) 

title ( 'Seeker head Yaw angle rate vs Time') 
meta jerryl, pause, clg 

plot (time (1: i-1) , betap(2 , (1: i-1) ) ) , grid, xlabel ( ' time' ) ,ylabe 
1 ( 'radians' ) 

title ( 'Seeker head Pitch angle rate vs Time') 
meta jerryl, pause, clg 

% commanded acceleration 

plot (time ( 1 : i) , acorn) , grid, title ( ' acceleration commanded' ) 
meta jerryl, pause, clg 

END 
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B. MISSILE-TARGET PROGRAM WITH TARGET OBSERVER 



clear 

! erase msltgt.met 
! erase jerryl.met 
clg 
% 

% 

% 

% This program simulates a 3 dimensional target/missile 

% engagement using basic proportional navigation with 

% an Luenberger observer and a ground target observer. 



% 

% 

% 

% 

% 

% 

% 

% 

% 

% 

% 



DEFINE STATES 
Missile 



xm 


missile 


xdm 


missile 


ym 


missile 


ydm 


- missile 


zm 


- missile 


zdm 


missile 



position in x direction 
velocity in x direction 
position in y direction 
velocity in y direction 
position in z direction 
velocity in z direction] 



Am= [ 0 100 
0 0 0 0 
0 0 0 1 
0 0 0 0 
0 0 0 0 
0 0 0 0 



0 0 
0 0 
0 0 
0 0 
0 1 
0 0 ]; 



Bm=[ 0 0 0 
10 0 
0 0 0 
0 10 
0 0 0 

0 0 1 ]; 

% 

% Seeker Head 

% 

% betap-[beta pitch - seeker head pitch angle 

% betad pitch - seeker head pitch angle rate 

% betay»[beta yaw - seeker head yaw angle 

% betad yaw - seeker head yaw angle rate] 

% 

% 



Asp=[0 1 ; -100 -20.0]; 
Asy= [ 0 1 ; -100 -20.0]; 
Bsp=[0 ; 100] ; 
Bsy=[0;100] ; 
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% 

% Autopilot 

% 

% gammamp= [ gammam pitch - pitch body angle 

% gammadm pitch - pitch body angle rate ] 

% gammamy= [ gammam yaw - yaw body angle 

% gammadm yaw - yaw body angle rate] 

% 

Ap=[-3 0 ; 0 -3]; 



Bp= [ 3 0 ; 0 3]; 

% Target 

% 



% 


ts=[xt 


- target 


position 


in 


X 


direction 


% 


xdt 


target 


velocity 


in 


X 


direction 


% 


yt 


- target 


position 


in 


y 


direction 


% 


ydt 


target 


velocity 


in 


y 


direction 


% 


zt 


target 


position 


in 


z 


direction 


% 


zdt 


- target 


velocity 


in 


z 


direction] 



% 

At= [ 0 10000 
0 0 0 0 0 0 
0 0 0 1 0 0 
0 0 0 0 0 0 
0 0 0 0 0 1 
000000 ]; 

Bt= [ 0 0 0 
10 0 
0 0 0 
0 10 
0 0 0 
ooi]; 

% 

% DISCRETE REPRESENTATION 

% 

dt=. 05; 

[phisp,delsp]= c2d(Asp,Bsp,dt) ; 
[phisy,delsy]= c2d(Asy, Bsy ,dt) ; 
[phim,delm]= c2d(Am,Bm,dt) ; 
[phia,dela]= c2d(Ap,Bp,dt) ; 
[phit,delt]= c2d(At,Bt,dt) ; 

tfinal= 20.0; 
kmax=tf inal/dt +1; 

% INITIAL VALUES 

% 
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% navigation ratio 

% 

NR= [4.0 0 ; 0 4.0]; 



% initial seeker head angles and angle rates 

% 

beta_pitchO =0.0; 
betad_pitch0=0 . 0 ; 
beta_yaw0 =0.0; 
betad_yawO =0.0; 

betap ( : , 1 ) = ( beta_pitchO 

betad_pitchO ] ; 
betay ( : , 1 ) = [ beta_yaw0 

betad_yawO ] ; 

% initial line-of-sight angular rates 

sigmad_pitch0=0 . 0 ; 
sigmad_yawO =0.0; 

sigmad( : , 1) =[sigmad_pitchO 
sigmad_yawO ] ; 

% initial observer angles and angle rates 

betae_pitch0=0 . 0 ; 
betade_pitch0=0 . 0 ; 
betae_yaw0=0 . 0 ; 
betaed_yaw0=0 . 0 ; 

% 

% 

betaep ( : , 1 ) = [ betaejpitchO 

betade_pitch0] ; 
betaey ( : , 1 ) = [ betae_yaw0 

betaed_yaw0 ] ; 

% 

% 

% initial missile body angle rates 

% 

gammadm_pitch0=0 . 0 ; 
gammadm_yawO =0.0; 

gammadm( : , 1) =[gammadm_pitch0 

gammadm yawO ] ; 

% 

% initial missile states 

% 

xmO =000.0; 
xdm0=2000 . 00; 
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ymO =000.0; 
ydm0=00; 
zmO =00; 
zdm0=00 ; 
ms ( : , 1) =(xmO 
xdmO 
ymO 
ydmO 
zmO 

zdmO ] ; 

% 

% initial target states 

% 

xtO =20000; 
xdt0=0000; 
ytO =050.0; 
ydt0=1000 . 0 ; 
ztO =050.0; 
zdt0=000. 0; 

ts ( : , 1) =[xtO 
xdtO 
ytO 
ydto 
ztO 

zdtO] ; 

% 

% initial range information 

% 

rxO=(xtO-xmO) ; 
rx ( 1) =rxO ; 
ryO= (ytO-ymO) ; 
ry(l)=ryO; 
rzO=(ztO-zmO) ; 
rz (1) =rzO; 

rmO=sqrt (xmO*2 + ym0*2 + zm0*2) ; 
rm ( 1 ) =rmO ; 

rtO=sqrt (xt0*2 + ytO*2 + ztO*2) ; 
rt (1) =rtO; 

rO =sqrt ( (xtO-xmO) *2 + (ytO-ymO) "2 + (ztO-zmO) "2) ; 
r ( 1) =rO ; 

% Initial Target Tracker values 

y ( ♦ / l)=[0;0;0] ; 

xhatl =[0;0;0;0;0;0] ; 

pkkml=[ 10*9 00000 
0 10*9 0000 
0 0 10*9 000 
000 10*9 0 0 
0000 10*9 0 
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0 0 0 0 0 10 A 9 ] ; 

rand( 'normal' ) ; 

h=[l 0 0 0 0 0 
0 0 1 0 0 0 
000010]; 

aT=[0 10000 
0 0 0 0 0 0 
0 0 0 1 0 0 
0 0 0 0 0 0 
0 0 0 0 0 1 
000000 ]; 

bT= [ 0 0 0 
10 0 
0 0 0 
0 10 
0 0 0 
0 0 1]; 

% 

dist=[0. 0013; 0.0513] ; 



% Continuous to Discrete 
% 

[phi,del]=c2d(aT,bT,dt) ; 

% 

% initial time 

% 

time (1) =0. 0; 

% 

% SIMULATION 

j=0; 

for (i=l:kmax-l) 

j-j+i; 

% missile and target flight path angles 

gammamjpitch(i)= s atan2 (ms (6, i) , sqrt (ms (2 , i) A 2+ms(4 , i) A 2) ) 
gammam_yaw(i) =atan2 (ms (4 , i) ,ms (2 , i) ) ; 

gammat_pitch(i) =atan2 (ts (6, i) , sqrt (ts(2 , i) A 2+ts(4, i) A 2) ) 
gamma t_yaw ( i ) =atan2 (ts(4,i) ,ts(2,i) ) ; 



% target velocity (magnitude) 

vt (i)=sqrt(ts(2,i) A 2+ts(4,i) A 2+ts(6,i) A 2) ; 
vm ( i) =sqrt (ms (2 , i) A 2+ms (4 , i) A 2+ms(6, i) A 2) ; 

% 

% line-of-sight angles 
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% 



sigma_pitch( i) =atan2 ( (ts (5, i) -ms (5, i) ) , sqrt ( (ts ( 1, i) -ms ( 1 , i) 
) * 2 . 

+ ( ts ( 3 , i ) -ms ( 3 , i ) ) A 2 ) ) ; 



sigma_yaw ( i) =atan2 ( (ts (3 , i) -ms (3 , i) ) , ( (ts ( l , i) -ms ( l , i) ) ) ) ; 

sigma ( : , i)=[sigma_pitch(i) 

sigma_yaw(i) ]; 



% update seeker head 

wk ( i ) =rand* . 9 ; 

% 

betap( : , i + l)=phisp*betap( : , i) 
delsp*sigma_pitch ( i) +dist*wk ( i) ; 

betay( : , i + l)=phisy*betay( : , i) 
delsy*sigma_yaw(i)+dist*wk(i) ; 

% 

% set up observer estimates 

% 

y0(l)=0.0; 
yl (1) =0. 0; 

% 

% 

gain=[-0 . 0031 ; 0 . 0299 ] ; 
c=(0 l]; 
f=phisp-gain*c; 
yO(i+l)=c*betap( : , i) ; 

betaep ( : , i+1) =f *betaep ( : , i) + gain*y0(i) 
delsp*sigma_pitch(i) ; 

yl(i+l)=c*betay( : , i) ; 
f l=phisy-gain*c; 

betaey ( : , i+l) =f l*betaey ( : , i) + gain*yl(i) 
delsy*sigma_yaw( i) ; 

% 

% 

% 

% seeker head angle rate vector 

% with uplinks at one second intervals 



+ 



+ 



+ 



+ 



if ( j==40) 

betad( : , i)=[ (betaep (2, i)+sigmad(l, i) ) 

(betaey(2, i)+sigmad(2, i) ) ] ; 
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j=0; 

else 

betad ( : , i ) = [ betaep ( 2 , i ) 

betaey ( 2 , i) ] ; 

end 



update autopilot 

acceleration commanded = vm(i) * gammadm(i) 

gammadm(i) = Nav constant * line of sight rate 

gammadm( : , i+1) =phia*gammadm( : , i) + dela*NR*betad( : , i) ; 

vm_pitch( i) =vm ( i) *cos (gammam yaw( i) -sigma_yaw( i) ) ; 
vm_yaw(i)=vm(i) *cos (gammam_pitch ( i) ) ; 



limit commanded acceleration to approximately 20 g's 

if (vm_yaw(i) *gammadm(2 , i) ) <=640. 0 

amc_y(i)=( v m_yaw(i) *gammadm(2, i) ) ; 

else 

amc_y ( i ) = 6 4 0 . 0 ; 

end 



if (vm_pitch(i) *gammadm( 1, i) ) <=640.0 

amc_p(i)=( v m_pitch(i) *gammadm(l, i) ) ; 

else 

amc_p ( i ) =64 0 . 0 ; 

end 

magnitude of commanded acceleration 
acom(i)=sqrt(amc_y(i) A 2 + amc_p(i) *2) ; 
missile acceleration vector components 



xddm(i) =(gammadm(l / i) *ms(6, i) *cos(gammam_yaw(i) ) ) ; 
yddm(i)=amc_y(i) *cos (gammam_yaw(i) ) ; 
zddm(i)= 2 unc_p(i) *cos(gammam_pitch(i) ) ; 

total missile acceleration magnitude 

am(i)=sqrt(xddm(i) A 2 + zddm(i) A 2 + yddm(i)*2 ); 

missile input vector 



<*> <*><*> X <#> <#> <#> <#> <#><*><#><#> <#><#>*> <*><#><#> <#><#> <*><*><#> <*><#><*> 



um= [ xddm ( i ) 
yddm ( i ) 
zddm(i) ] ; 

update missile 

ms ( : , i+1) =phim*ms ( : , i) + delm*um; 
set target acceleration components 
if (r (i) <=00000.0) 

xddt(i)“ -5.0*32.2*sin(sigma_yaw(i) ) ; 
yddt(i)= 3 . 5*32 . 2*cos (sigma_yaw(i) ) ; 
zddt ( i) = 3 . 5*32 . 2*cos (sigma_pitch ( i) ) ; 

else 

xddt ( i ) =0 . 0 ; 
yddt(i)=0.0; 
zddt ( i) =0 . 0; 

end 

target acceleration magnitude 

at ( i) =sqrt (xddt ( i) *2 + yddt(i)*2 + zddt(i) A 2); 

target input vector 

ut=[xddt(i) 
yddt (i) 
zddt(i) ] ; 

update target 

ts ( : , i+1) =phit*ts ( : , i) + delt*ut; 



missile and target trajectory data 

missile(i, : ) “[ms (1, i) ms(3,i) ms(5,i)]; 
target (i, : ) »[ts(l,i) ts(3,i) ts(5,i)]; 

Target traker 

Input initial conditions 



(: ,i)-[ts(l,i) ;ts(2, i) ;ts(3,i) ;ts(4,i) ;ts(5,i) ;ts(6,i) ] ; 
Plant noise covariance 

Ql“ [ 160'“ 2 0 0; 0 160 A 2 0;0 0 160 A 2]; 
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<*><#><#> 



% Measurement noise covariance 



R= [500*2 0 0 
0 125*2 0 
0 0 12 * 2 ] ; 

% Input equations 

Q=[ 160 *2; 160* 2; 160 *2] ; 

nu=[sqrt (500) *2 ; sqrt ( 125*2 ) ;sqrt(12*2) ]*rand; 
w(:,i) =sqrt (Q) *rand; 



% Update the plant states 
% 

x( : , i+l)=phi*x( : , i) +del*w( : , i) ; 
y ( : , i+l)=h*x( : , i+1) +nu; 

Compute Kalman gains 

G = pkkml*h' *inv(h*pkkml*h'+R) ; 

% 

% Define gain matrix for plotting 

g( : / i) = [ G ( 1 : 2 , 1 ) ; G ( 3 : 4 , 2 ) ;G(5:6,3) ] ; 

% 

% Measurement step 

xhat (:,i+l) =xhatl+(G * (y ( : , i+1) -h*xhatl) ) ; 
pkk =(eye(6)-G *h)*pkkml; 

% 

% Movement step 

xhat l=phi* xhat ( : , i+1) ; 
pkkml=phi*pkk*phi'+ del*Ql*del'; 

% compute the range deviation 

xhatr =[xhat(l, i) ;xhat(3, i) ;xhat(5, i) ] ; 
yr =y(:,i); 

xr(:,i)» (yr - xhatr); 

% compute the velocity deviation 

xhatv =[xhat (2 , i) ;xhat (4 , i) ;xhat (6, i) ] ; 
yv =[ts(2,i) ;ts(4,i) ;ts(6,i) ] ; 

xv(:,i)= (yv - xhatv); 

% 
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% Transformation to body axes 



all=cos (gammam_pitch ( i) ) *cos (gammam_yaw( i) ) ; 
al2=cos (gammam_pitch ( i) ) *sin(gammam_yaw(i) ) ; 
al3=(-sin(gammam_pitch( i) ) ) ; 
a21=(-sin(gammam_yaw(i) ) ) ; 
a22=cos (gammam_yaw( i) ) ; 
a23=0 . 0 ; 

a31=sin(gammam_pitch(i) ) *cos (gammam_yaw ( i) ) ; 
a32=sin(gammam_pitch(i) ) *sin(gammam_yaw( i) ) ; 
a33=cos (gammam_pitch ( i) ) ; 

% velocity transformation 

vbx(i)=xv(l, i) *all+xv(2, i) *al2+xv(3, i) *al3 ; 
vby (i) =xv( 1, i) *a21+xv(2 , i) *a22+xv(3 , i) *a23 ; 
vbz (i) =xv( 1, i) *a3 1+xv ( 2 , i) *a32+xv(3 , i) *a33 ; 

% range transformation 

rbx( i) =xr ( 1, i) *all+xr (2 , i) *al2+xr (3 , i) *al3 ; 
rby (i) =xr (1, i) *a21+xr (2 , i) *a22+xr (3 , i) *a23 ; 
rbz (i) =xr (1, i) *a31+xr (2 , i) *a32+xr (3 , i) *a33 ; 

% line-of-sight transformation 

sll=cos (sigma_pitch(i) ) *cos (sigma_yaw( i) ) ; 
sl2=cos (sigma_pitch ( i) ) *sin(sigma_yaw(i) ) ; 
sl3=(-sin(sigma_pitch(i) ) ) ; 
s21=(-sin(sigma_yaw(i) ) ) ; 
s22=cos (sigma_yaw( i) ) ; 
s23=0 . 0 ; 

s31=sin(sigma_pitch(i) ) *cos (sigma_yaw( i) ) ; 
s32=sin(sigma_pitch(i) ) *sin(sigma_yaw(i) ) ; 
s33=cos(sigma_pitch(i) ) ; 

% velocity transformation 

vsx(i)=vbx(i) *sll+vby (i) *sl2+vbz (i) *sl3; 
vsy(i)=vbx(i) *s21+vby(i) *s22+vbz(i) *s23; 
vsz ( i) »vbx ( i) *s31+vby(i) *s32+vbz(i) *s33; 

rsx(i)=rbx(i) *sll+rby(i) *sl2+rbz (i) *sl3 ; 
rsy (i)=rbx(i) *s21+rby (i) *s22+rbz (i) *s23; 
rsz (i)=rbx(i) *s31+rby(i) *s32+rbz (i) *s33; 

% compute the magnitude of the range deviation 

rd(i)=sqrt(rsx(i) ^2+rsy (i) *2+rsz(i) A 2) ; 

% generate the angular rate deviations 
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sigmad_pitch(i+l)= vsy(i) / (r (i) +rd(i) ) ; 
sigmad_yaw(i+l) = vsz (i) / (r(i) +rd(i) ) ; 

sigmad( : , i+1) =[sigmad_pitch(i) 

s igmad_yaw ( i ) ] ; 



% update range information 
% 



r (i + 1) =sqrt ( (ts (1, i + 1) -is(l , i + l) ) ~ 2 + 

(ts (3, i+1) -ms (3, i+1) ) ~2. . . 

+ (ts(5, i+1) -ms(5, i+1) ) A 2) ; 
rm(i+l) =sqrt (ms ( 1, i+1 ) "2 + ms (3, i+1) A 2 + ms (5, i+1) "2) ; 
rt (i+1) =sqrt (ts (1 , i+1) A 2 + ts (3 , i+1) A 2 + ts (5, i+1) "2) ; 

% update time vector 

% 

time ( i+1) =time(i) + dt; 

% 

% check to see if engagement is at closest point of 

approach 
% 

if (r(i)<r(i+l) ) , break, end 

end; 

% 

% PLOT RESULTS 

% 

% range data 

% 

plot(time,r) , grid, xlabel ( 'TIME' ) , y label (' FEET' ) 
title ('RANGE VS TIME') 

text( 7,20000, ( 'rng=' , num2str (r ( i) ), ' feet']); 
meta msltgt, pause, clg 

subplot(211) ,plot (time,rm) , grid, xlabel( 'TIME' ) , y label ( 'FEET' ) 
title ('MISSILE RANGE VS TIME'); 

subplot(212) , plot (time, rt) , grid, xlabel( 'TIME' ) , y label ( 'FEET' ) 
title ('TARGET RANGE VS TIME'); 
meta msltgt, pause, clg 
% 

% 

% missile and target velocity 

% 

%plot (time(l: i) ,vm) , grid, xlabel ( 'TIME' ) , y label ( 'FEET/ SEC' ) 
%title( 'MISSILE VELOCITY VS TIME'); 

%meta msltgt, pause, clg 

%plot (time (1 : i) ,vt) , grid, xlabel ( 'TIME' ) , y label ( 'FEET/ SEC' ) 
%title( 'TARGET VELOCITY VS TIME'); 

Imeta msltgt, pause, clg 
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% 



missile and target acceleration 



plot (time (1: i) , am) , grid, xlabel ( 'TIME' ) , y label ( ' FEET/ SEC* 2 ' ) 
title ( 'MISSILE ACCELERATION vs TIME') 
meta msltgt, pause, clg 

plot (time(l: i) , at) , grid, xlabel ( 'time' ) ,y label ( ' feet/ sec'' 2 ' ) 
title ( 'TARGET ACCELERATION vs TIME') 
meta msltgt, pause, clg 

% seeker head angle rates 

plot (time (1 : i— 1 ) , betaey (2 , ( 1: i-1) ) ) , grid, xlabel ( ' time' ) , ylab 
el ( 'radians' ) 

title ( 'SEEKER HEAD YAW ANGLE RATE vs TIME') 
meta msltgt, pause, clg 

plot (time(l: i-1) , betaep(2 , (1: i-1) ) ) , grid, xlabel ( ' time' ) ,ylab 
el ( 'radians' ) 

title ('SEEKER HEAD PITCH ANGLE RATE vs TIME') 
meta msltgt , pause , clg 

% commanded acceleration 

plot(time(l: i) , acorn) , grid, title ( 'ACCELERATION COMMANDED') 
xlabel ( 'TIME' ) ,ylabel ( ' FEET/SEC~2 ' ) 
meta msltgt , pause, clg 



% Kalman estimates 

plot (time(l: i-l) ,x(l, (1: i-1) ) , ' , time (1: i-1) , xhat (1, ( 1: i-l) 

),'+') 

xlabel ( 'TIME') , ylabel ( ' FEET' ) ,grid 

title ( 'TARGET POSITION in the X DIRECTION vs TIME') 
gtext ( ' ** ACTUAL TARGET STATE ') 
gtext ( ' ++ KALMAN FILTER ESTIMATE') 
meta msltgt , pause, clg 

plot (time (1: i-l) ,x(2, (1: i-1) ) , '*' , time (1: i-1) , xhat (2 , (1: i-1) 

) ,'+') 

xlabel ('TIME') , ylabel ( 'FEET/SEC' ) ,grid 

title ( 'TARGET VELOCITY in the X DIRECTION vs TIME') 

gtext ('** ACTUAL TARGET STATE ') 

gtext ('++ KALMAN FILTER ESTIMATE') 

meta msltgt, pause, clg 

plot ( time (l:i-50) ,g(l, (l:i-50) ) , '*' ,time(l: i-50) ,g(2, (l:i-50 
) ) , '+') 

title ('KALMAN GAINS in the X DIRECTION' ), grid, xlabel ( 'TIME' ) , 
gtext ('*** (POSITION GAIN) ' ) 
gtext ( ' +++ (VELOCITY GAIN) ' ) 
meta msltgt, pause, clg 
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% 



Kalman estimates 



plot (time(l: i-1) , x(3, ( 1 : i— 1 ) ) , '*' , time(l:i-i) ,xhat(3, (l:i-i) 

) , '+') 

xlabel ( 'TIME' ) , ylabel ( 'FEET' ) ,grid 

title ( 'TARGET POSITION in the Y DIRECTION vs TIME') 
gtext ( ' ** ACTUAL TARGET STATE ') 
gtext ( ' ++ KALMAN FILTER ESTIMATE') 
meta msltgt , pause, clg 

plot (time (1: i-1) ,x(4, (1: i-1) ) , '*' , time(l: i-1) ,xhat (4, (1: i-1) 
) , '+') 

xlabel ('TIME') , ylabel ( ' FEET/SEC' ) ,grid 

title( 'TARGET VELOCITY in the Y DIRECTION VS TIME') 

gtext ('** ACTUAL TARGET STATE ') 

gtext ('++ KALMAN FILTER ESTIMATE') 

meta msltgt, pause, clg 

plot (time (1 : i-50) ,g(3, (1: i-50) ) , ' * ' , time ( 1 : i-50) ,g(4, (l:i-50 

)),'+') 

title ( 'KALMAN GAINS in the Y DIRECTION' ), grid, Xlabel ( 'TIME' ) , 
gtext ('*** (POSITION GAIN) ' ) 
gtext ( ' +++ (VELOCITY GAIN) ' ) 
meta msltgt , pause, clg 

plot (time (1 : i-1) ,x(5, ( 1: i-1) ) , '*' , time ( 1: i-1) , xhat (5 , ( 1 : i-1) 

),'+') 

xlabel ( 'TIME' ) , ylabel (' FEET' ) ,grid 

title( 'TARGET POSITION in the Z DIRECTION vs TIME') 
gtext ('** ACTUAL TARGET STATE ') 
gtext ('++ KALMAN FILTER ESTIMATE') 
meta msltgt , pause, clg 

plot (time (1 : i-1) ,x(6, (1: i-1) ) time ( 1: i-1) , xhat (6, (1 : i-1) 

),'+') 

xlabel ('TIME') , ylabel ( 'FEET/SEC') ,grid 

title( 'TARGET VELOCITY in the Z DIRECTION vs TIME') 

gtext ('** ACTUAL TARGET STATE ') 

gtext ('++ KALMAN FILTER ESTIMATE') 

meta msltgt, pause, clg 

plot (time(l: i-50) ,g(5, (1: i-50) ) , '*', time (l: i-50) ,g(6, (1: i-50 
) ) , '+') 

title ( 'KALMAN GAINS in the Z DIRECTION' ), grid, xlabel ( 'TIME' ) , 

gtext ( ' *** (POSITION GAIN) ' ) 

gtext ( ' +++ (VELOCITY GAIN) ' ) 

meta msltgtl , pause, clg 

end 
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